function BootAll

clear;clc;
global RMarket_M RMarket_Q RMarket_ADec RMarket_AJun RFree_M RFree_Q RFree_ADec RFree_AJun MktVar_M MktVar_Q MktVar_ADec MktVar_AJun yyyymm yyyyq yyyy SimpleLog BootSample BootISstruct BootOOSstruct rootdir;

rootdir = '...';

load(strcat(rootdir,'Data\P00Data'));
load(strcat(rootdir,'Data\Market'));

BootISstruct = struct('beta',[], 'tOls',[], 'tNwest',[], 'adjr2',[]);
BootOOSstruct = struct('R2',[], 'DMAE',[], 'MSET',[], 'MSEF',[], 'ENCT',[], 'ENCF',[], 'ENClambd',[], ...
    'R2_CT',[], 'DMAE_CT',[], 'MSET_CT',[], 'MSEF_CT',[], 'ENCT_CT',[], 'ENCF_CT',[], 'ENClambd_CT',[], ...
    'DiffUtil1',[], 'DiffUtil1t',[], 'DiffUtil2',[], 'DiffUtil2t',[], 'DLSt', [], 'DLSSRt', []);

NBoot = 10000;
gam = 5;
NVARS = 46;

% vari = [1 9 25 28] (pce, ogap, tchi, shtint) are regular reconstruction
% vari = [13 17] (sntm, fbm) are PLS
% vari = 22 is rsvix
% vari = 46 is cay

% DGP is NOT AR1 for vari = [1 9 23 24 25 46] (pce, ogap, gpce, gip, tchi, cay)

% In all, special cases
%  1 (pce)      Regular reconstruction / Not AR1
%  9 (ogap)     Regular reconstruction / Not AR1
% 13 (sntm)     PLS
% 17 (fbm)      PLS
% 22 (rsvix)    ??
% 23 (gpce)     Not AR1
% 24 (gip)      Not AR1
% 25 (tchi)     Regular reconstruction / Not AR1
% 28 (shtint)   Regular reconstruction
% 46 (cay)      Not AR1

% Examples (not special cases)
% vari = 10; %monthly
% vari = 8; %quarterly
% vari = 26; %annual

for BootSample = [1 3]
    if BootSample==1
        load(strcat(rootdir,'Analysis\Sample1_Full\ResultsFile'));
    elseif BootSample==3
        load(strcat(rootdir,'Analysis\Sample3_BegAuthor_EndFull\ResultsFile'));
    end

    for SimpleLog = 1:2
        fprintf('BootSample = %d, SimpleLog = %d\n', BootSample,SimpleLog);
        AllBootStatsIS_M = repmat(BootISstruct,NBoot,NVARS); AllBootStatsIS_Q  = repmat(BootISstruct,NBoot,NVARS); AllBootStatsIS_ADec = repmat(BootISstruct,NBoot,NVARS); AllBootStatsIS_AJun = repmat(BootISstruct,NBoot,NVARS);
        AllBootStatsOOS_M = repmat(BootOOSstruct,NBoot,NVARS); AllBootStatsOOS_Q = repmat(BootOOSstruct,NBoot,NVARS); AllBootStatsOOS_ADec = repmat(BootOOSstruct,NBoot,NVARS); AllBootStatsOOS_AJun = repmat(BootOOSstruct,NBoot,NVARS);

        if SimpleLog==1
            RMarket_M = CRSPSP_M; RFree_M = Rf_M;  MktVar_M = MktStd(:,1).^2;
            RMarket_Q = CRSPSP_Q(3:3:end); RFree_Q = Rf_Q(3:3:end); MktVar_Q = 3*(MktStd(:,1).^2); MktVar_Q = MktVar_Q(3:3:end);
            RMarket_ADec = CRSPSP_A(12:12:end); RFree_ADec = Rf_A(12:12:end); MktVar_ADec = 12*(MktStd(:,1).^2); MktVar_ADec = MktVar_ADec(12:12:end);
            RMarket_AJun = CRSPSP_A(6:12:end); RFree_AJun = Rf_A(6:12:end); MktVar_AJun = 12*(MktStd(:,1).^2); MktVar_AJun = MktVar_AJun(6:12:end);
        else
            RMarket_M = log(1+CRSPSP_M); RFree_M = log(1+Rf_M); MktVar_M = MktStd(:,2).^2;
            RMarket_Q = log(1+CRSPSP_Q(3:3:end)); RFree_Q = log(1+Rf_Q(3:3:end)); MktVar_Q = 3*(MktStd(:,2).^2); MktVar_Q = MktVar_Q(3:3:end);
            RMarket_ADec = log(1+CRSPSP_A(12:12:end)); RFree_ADec = log(1+Rf_Q(12:12:end)); MktVar_ADec = 12*(MktStd(:,2).^2); MktVar_ADec = MktVar_ADec(12:12:end);
            RMarket_AJun = log(1+CRSPSP_A(6:12:end)); RFree_AJun = log(1+Rf_Q(6:12:end)); MktVar_AJun = 12*(MktStd(:,2).^2); MktVar_AJun = MktVar_AJun(6:12:end);
        end

        if BootSample==1
            varRange = 1:NVARS;
        else
            varRange = 1:29;
        end

        for vari = varRange
            disp(vari);

            if vari==1 %Regular reconstruction of "pce" / DGP is not AR1
                BootStatsIS_M = repmat(BootISstruct,NBoot,1); BootStatsOOS_M = repmat(BootOOSstruct,NBoot,1);

                [TBeg, TEnd, OOSBeg, MeanVar, MeanRf] = CreateParams(vari,'Quarterly');
                RM = RMarket_Q; RF = RFree_Q;
                X = readmatrix(strcat(rootdir,'Data\P01_Data.xlsx'), 'Sheet','Quarterly', 'Range','a3:e302');
                yq = X(:,1); PCE = log(X(:,3) + X(:,4)); % these data are from 1947:1 to 2021:4
                T1 = find(yyyyq==yq(1)); T2 = TEnd(1);
                Rm = RM(T1+1:T2); Rf = RF(T1+1:T2);
                tmp = find(yq==yyyyq(TEnd(1)));
                PCE(tmp+1:end) = []; % these data are from 1947:1 to 2021:4 (or TEnd)
                gPCE = PCE(2:end) - PCE(1:end-1); % these data are from 1947:2 to 2021:4 (or TEnd)
                clear T1 T2 yq tmp;

                [gBoot_Q, RmBoot_Q, RfBoot_Q] = BootGenData(gPCE, Rm, Rf, NBoot); %Data from 1947:2 to 2021:4 (or TEnd)
                PCEBoot_Q = PCE(1) + cumsum(gBoot_Q); %Data from 1947:2 to 2021:4 (or TEnd)
                % I need data from 1947:1 to 2021:4 (or TEnd); Take the first observation to be always actual for PCE
                PCEBoot_Q = [repmat(PCE(1), 1, NBoot); PCEBoot_Q];
                RmBoot_Q = [repmat(mean(Rm), 1, NBoot); RmBoot_Q];
                RfBoot_Q = [repmat(mean(Rf), 1, NBoot); RfBoot_Q];

                T = size(RmBoot_Q, 1); K = 24;
                XBootOOS_Q = NaN(T, T, NBoot);
                for b = 1:NBoot
                    for t = 40:T
                        PCE_filter = NaN(T, 1);
                        Y = PCEBoot_Q(K+4:t, b); X = [ones(length(Y),1) PCEBoot_Q(4:t-K, b) PCEBoot_Q(3:t-K-1, b) PCEBoot_Q(2:t-K-2, b) PCEBoot_Q(1:t-K-3, b)];
                        bb = X\Y; PCE_filter(K+4:t) = Y - X*bb;
                        XBootOOS_Q(:, t, b) = PCE_filter;
                    end
                end
                XBootIS_Q = squeeze(XBootOOS_Q(:, end, :));
                % remove the first K+3 = 27 observations so that the data correspond to 1953:4 to 2021:4 (or TEnd) sample
                RmBoot_Q(1:27, :) = []; RfBoot_Q(1:27, :) = []; XBootIS_Q(1:27, :) = []; XBootOOS_Q = XBootOOS_Q(28:end, 28:end, :);

                [BootStatsIS_Q, BootStatsOOS_Q, BootStatsIS_ADec, BootStatsOOS_ADec, BootStatsIS_AJun, BootStatsOOS_AJun] = ...
                    RunQuarterly(XBootIS_Q, XBootOOS_Q, RmBoot_Q, RfBoot_Q, OOSBeg, gam, MeanVar, MeanRf, NBoot);

            elseif vari==9 %Regular reconstruction of "ogap"  / DGP is not AR1
                [TBeg, TEnd, OOSBeg, MeanVar, MeanRf] = CreateParams(vari,'Monthly');
                RM = RMarket_M; RF = RFree_M;
                Rm = RM(TBeg(1):TEnd(1)); Rf = RF(TBeg(1):TEnd(1));
                IP = readmatrix(strcat(rootdir,'Data\P01_Data.xlsx'), 'Sheet','Monthly', 'Range','c86:c1238'); % read from 1925:12 to 2021:12
                % convert to data from TBeg to TEnd
                IP = IP(TBeg:TEnd+1);
                IP = log(IP);
                gIP = IP(2:end) - IP(1:end-1); % from 1926:01 to 2021:12 (TBeg to TEnd)
                [gBoot_M, RmBoot_M, RfBoot_M] = BootGenData(gIP, Rm, Rf, NBoot); % bootstrap growthIP
                IPBoot_M = IP(1) + cumsum(gBoot_M); % convert to IP  ; from 1926:01 to 2021:12
                T = size(IPBoot_M,1);

                XRegr = [ones(1,T); (1:T); (1:T).^2]';
                XBootOOS_M = NaN(T, T, NBoot);
                for b = 1:NBoot
                    for t = 60:T
                        y = IPBoot_M(1:t, b); x = XRegr(1:t, :);
                        XBootOOS_M(1:t, t, b) = y - x*(x\y);
                    end
                end
                XBootIS_M = squeeze(XBootOOS_M(:, end, :));
                [BootStatsIS_M, BootStatsOOS_M, BootStatsIS_Q, BootStatsOOS_Q, BootStatsIS_ADec, BootStatsOOS_ADec, BootStatsIS_AJun, BootStatsOOS_AJun] = ...
                    RunMonthly(XBootIS_M, XBootOOS_M, RmBoot_M, RfBoot_M, OOSBeg, gam, MeanVar, MeanRf, NBoot);

            elseif vari==25 % Paper #22 Regular reconstruction of "tchi"  / DGP is not AR1
                [TBeg, TEnd, OOSBeg, MeanVar, MeanRf] = CreateParams(vari,'Monthly');
                RM = RMarket_M; RF = RFree_M;
                Rm = RM(TBeg(1):TEnd(1)); Rf = RF(TBeg(1):TEnd(1));

                % Have to generate bootstrap data here as this is a special case
                % Basically draw 1s/0s from XALL randomly; then do PCA
                load(strcat(rootdir,'Data\P22_Data'), 'XALL', 'ym'); XALL(1:12,:) = []; ym(1:12) = [];
                tmp = find(ym==yyyymm(TEnd(1)));
                XALL(tmp+1:end,:) = [];
                clear ym tmp;
                T = size(XALL,1);
                muRm = mean(Rm); eRm = Rm - muRm;
                muRf = mean(Rf); eRf = Rf - muRf;
                RmBoot_M = NaN(T, NBoot); RfBoot_M = NaN(T, NBoot); XBootOOS_M = NaN(T, T, NBoot);
                TIndex = randi([1 T], T, NBoot);
                for b = 1:NBoot
                    Tb = TIndex(:,b);
                    RmBoot_M(:,b) = muRm + eRm(Tb);
                    RfBoot_M(:,b) = muRf + eRf(Tb);
                    Xb = XALL(Tb,:);
                    for t = 60:T
                        [~,TECH] = pca(Xb(1:t,:));
                        XBootOOS_M(1:t, t, b) = TECH(:,1);
                    end
                end
                XBootIS_M = squeeze(XBootOOS_M(:, end, :));
                [BootStatsIS_M, BootStatsOOS_M, BootStatsIS_Q, BootStatsOOS_Q, BootStatsIS_ADec, BootStatsOOS_ADec, BootStatsIS_AJun, BootStatsOOS_AJun] = ...
                    RunMonthly(XBootIS_M, XBootOOS_M, RmBoot_M, RfBoot_M, OOSBeg, gam, MeanVar, MeanRf, NBoot);

            elseif vari==28 % Paper #25 Regular reconstruction of "shtint"
                [TBeg, TEnd, OOSBeg, MeanVar, MeanRf] = CreateParams(vari,'Monthly');
                RM = RMarket_M; RF = RFree_M;
                Rm = RM(TBeg(1):TEnd(1)); Rf = RF(TBeg(1):TEnd(1));

                load(strcat(rootdir,'Data\P25_Data'), 'EWSI'); EWSI = EWSI(TBeg(1):TEnd(1));
                EWSI = log(EWSI); % valid data are from 1973:12
                gEWSI = EWSI(2:end) - EWSI(1:end-1); % these data are from 1974:01
                [gBoot_M, RmBoot_M, RfBoot_M] = BootGenDataAR1(gEWSI, Rm(2:end), Rf(2:end), NBoot); % bootstrap growthEWSI
                EWSIBoot_M = EWSI(1) + cumsum(gBoot_M); % convert to EWSI

                % I need data from 1973:12 ; Take the first observation to be always actual
                EWSIBoot_M = [repmat(EWSI(1), 1, NBoot); EWSIBoot_M];
                RmBoot_M = [repmat(mean(Rm), 1, NBoot); RmBoot_M];
                RfBoot_M = [repmat(mean(Rf), 1, NBoot); RfBoot_M];
                T = size(RmBoot_M,1); XRegr = [ones(1,T); (1:T)]';
                XBootOOS_M = NaN(T, T, NBoot);
                for b = 1:NBoot
                    for t = 60:T
                        yreg = EWSIBoot_M(1:t,b); xreg = XRegr(1:t,:);
                        sii = yreg - xreg*(xreg\yreg);
                        XBootOOS_M(1:t, t, b) = sii/std(sii);
                    end
                end
                XBootIS_M = squeeze(XBootOOS_M(:, end, :));
                [BootStatsIS_M, BootStatsOOS_M, BootStatsIS_Q, BootStatsOOS_Q, BootStatsIS_ADec, BootStatsOOS_ADec, BootStatsIS_AJun, BootStatsOOS_AJun] = ...
                    RunMonthly(XBootIS_M, XBootOOS_M, RmBoot_M, RfBoot_M, OOSBeg, gam, MeanVar, MeanRf, NBoot);

            elseif vari==46 %Regular reconstruction of "cay"
                BootStatsIS_M = repmat(BootISstruct,NBoot,1); BootStatsOOS_M = repmat(BootOOSstruct,NBoot,1);

                [TBeg, TEnd, OOSBeg, MeanVar, MeanRf] = CreateParams(vari,'Quarterly');
                RM = RMarket_Q; RF = RFree_Q;
                Rm = RM(TBeg(1):TEnd(1)); Rf = RF(TBeg(1):TEnd(1));
                x = load('../../../../Comprehensive/2021/Data/cay_q_21Q4_me.txt');
                C = x(:,2); A = x(:,3); Y = x(:,4); clear x; % valid data are from 1952:1 to 2021:4
                gCAY = [C(2:end) - C(1:end-1) A(2:end) - A(1:end-1) Y(2:end) - Y(1:end-1)];  % these data are from 1952:2 to 2021:4
                [gBoot_Q, RmBoot_Q, RfBoot_Q] = BootGenData(gCAY, Rm(2:end), Rf(2:end), NBoot); %data are from 1952:2 to 2021:4
                CBoot_Q = C(1) + cumsum(gBoot_Q(:,:,1)); ABoot_Q = A(1) + cumsum(gBoot_Q(:,:,2)); YBoot_Q = Y(1) + cumsum(gBoot_Q(:,:,3));
                % I need data from 1952:1 to 2021:4; Take the first observation to be always actual for PCE
                CBoot_Q = [repmat(C(1), 1, NBoot); CBoot_Q];
                ABoot_Q = [repmat(A(1), 1, NBoot); ABoot_Q];
                YBoot_Q = [repmat(Y(1), 1, NBoot); YBoot_Q];
                RmBoot_Q = [repmat(mean(Rm), 1, NBoot); RmBoot_Q];
                RfBoot_Q = [repmat(mean(Rf), 1, NBoot); RfBoot_Q];

                T = size(RmBoot_Q, 1);
                XBootOOS_Q = NaN(T, T, NBoot);
                for b = 1:NBoot
                    DA = [NaN; ABoot_Q(2:end,b) - ABoot_Q(1:end-1,b)];
                    DY = [NaN; YBoot_Q(2:end,b) - YBoot_Q(1:end-1,b)];
                    for t = 21:103
                        YVAR = CBoot_Q(4:t-2,b);
                        XVAR = [ones(length(YVAR),1) ABoot_Q(4:t-2,b) YBoot_Q(4:t-2,b)];
                        for k = -2:2
                            XVAR = [XVAR DA(4+k:t-2+k) DY(4+k:t-2+k)];
                        end
                        B = XVAR\YVAR;
                        cay_t = CBoot_Q(:,b) - B(1) - B(2)*ABoot_Q(:,b) - B(3)*YBoot_Q(:,b);
                        cay_t = cay_t - mean(cay_t);
                        XBootOOS_Q(1:t, t ,b) = cay_t(1:t);
                    end
                    for t = 104:T
                        YVAR = CBoot_Q(10:t-8,b);
                        XVAR = [ones(length(YVAR),1) ABoot_Q(10:t-8,b) YBoot_Q(10:t-8,b)];
                        for k = -8:8
                            XVAR = [XVAR DA(10+k:t-8+k) DY(10+k:t-8+k)];
                        end
                        B = XVAR\YVAR;
                        cay_t = CBoot_Q(:,b) - B(1) - B(2)*ABoot_Q(:,b) - B(3)*YBoot_Q(:,b);
                        cay_t = cay_t - mean(cay_t);
                        XBootOOS_Q(1:t, t ,b) = cay_t(1:t);
                    end
                end
                XBootIS_Q = squeeze(XBootOOS_Q(:, end, :));

                [BootStatsIS_Q, BootStatsOOS_Q, BootStatsIS_ADec, BootStatsOOS_ADec, BootStatsIS_AJun, BootStatsOOS_AJun] = ...
                    RunQuarterly(XBootIS_Q, XBootOOS_Q, RmBoot_Q, RfBoot_Q, OOSBeg, gam, MeanVar, MeanRf, NBoot);

            elseif vari==13 %Paper #12 "sntm" PLS
                [TBeg, TEnd, OOSBeg, MeanVar, MeanRf] = CreateParams(vari,'Monthly');
                RM = RMarket_M; RF = RFree_M;
                Rm = RM(TBeg(1):TEnd(1)); Rf = RF(TBeg(1):TEnd(1));

                load(strcat(rootdir,'Data\P12_Data'), 'BWVARS', 'yearmo');
                BWVARS = squeeze(BWVARS(:,end,:)); % there are actually different vintages. But for bootstrap use the latest vintage
                tmp = find(yearmo==yyyymm(TEnd(1)));
                BWVARS(tmp+1:end,:) = [];
                clear tmp;
                [XPrePLSBoot_M, RmBoot_M, RfBoot_M] = BootGenDataAR1(BWVARS, Rm, Rf, NBoot);
                [BootStatsIS_M, BootStatsOOS_M, BootStatsIS_Q, BootStatsOOS_Q, BootStatsIS_ADec, BootStatsOOS_ADec, BootStatsIS_AJun, BootStatsOOS_AJun] = ...
                    RunMonthlyPLS(XPrePLSBoot_M, RmBoot_M, RfBoot_M, OOSBeg, gam, MeanVar, MeanRf, NBoot);

            elseif vari==17 %Paper #16 "fbm" PLS
                [TBeg, TEnd, OOSBeg, MeanVar, MeanRf] = CreateParams(vari,'Monthly');
                RM = RMarket_M; RF = RFree_M;
                Rm = RM(TBeg(1):TEnd(1)); Rf = RF(TBeg(1):TEnd(1));

                load(strcat(rootdir,'Data\P16_Data'), 'BMAll');
                BMAll = BMAll(TBeg(1):TEnd(1),:); T = size(BMAll, 1);
                BMAll = BMAll./repmat(std(BMAll,'omitnan'),T,1);
                [XPrePLSBoot_M, RmBoot_M, RfBoot_M] = BootGenDataAR1(BMAll, Rm, Rf, NBoot);
                [BootStatsIS_M, BootStatsOOS_M, BootStatsIS_Q, BootStatsOOS_Q, BootStatsIS_ADec, BootStatsOOS_ADec, BootStatsIS_AJun, BootStatsOOS_AJun] = ...
                    RunMonthlyPLS(XPrePLSBoot_M, RmBoot_M, RfBoot_M, OOSBeg, gam, MeanVar, MeanRf, NBoot);

            elseif vari==22 %Paper #20 "rsvix
                [TBeg, TEnd, OOSBeg, MeanVar, MeanRf] = CreateParams(vari,'Monthly');
                RM = RMarket_M; RF = RFree_M;

                % run montly stuff
                XIS = VARS(vari).DataIS{1}; X = XIS(TBeg(1):TEnd(1)); Rm = RM(TBeg(1):TEnd(1)); Rf = RF(TBeg(1):TEnd(1));
                [XBootIS_M, RmBoot_M, RfBoot_M] = BootGenDataAR1(X, Rm, Rf, NBoot);
                T = size(XBootIS_M,1); XBootOOS_M = reshape(repmat(XBootIS_M, [T 1]), T, T, NBoot);
                YBoot_M = RmBoot_M - RfBoot_M;
                BootStatsIS_M = BootIS(XBootIS_M, YBoot_M);
                PredSgn = sign(extractfield(BootStatsIS_M,'beta'))';
                BootStatsOOS_M = BootOOS(XBootOOS_M, YBoot_M, OOSBeg(1), PredSgn, gam, MeanVar(1), MeanRf(1));

                % run quarterly stuff - Have to bootstrap data again as the X variable is different
                XIS = VARS(vari).DataIS{2}; X = XIS(TBeg(1):TEnd(1)); Rm = RM(TBeg(1):TEnd(1)); Rf = RF(TBeg(1):TEnd(1));
                [XBootIS_M, RmBoot_M, RfBoot_M] = BootGenDataAR1(X, Rm, Rf, NBoot);
                T = size(XBootIS_M,1); XBootOOS_M = reshape(repmat(XBootIS_M, [T 1]), T, T, NBoot);
                XBootIS_Q = XBootIS_M(3:3:end, :);
                XBootOOS_Q = XBootOOS_M(3:3:end,3:3:end,:);
                T = size(XBootIS_Q,1); RmBoot_Q = NaN(T, NBoot); RfBoot_Q = NaN(T, NBoot);
                for t = 1:T
                    RmBoot_Q(t, :) = prod(1+RmBoot_M(3*t-2:3*t, :)) - 1;
                    RfBoot_Q(t, :) = prod(1+RfBoot_M(3*t-2:3*t, :)) - 1;
                end
                YBoot_Q = RmBoot_Q - RfBoot_Q;
                BootStatsIS_Q = BootIS(XBootIS_Q, YBoot_Q);
                PredSgn = sign(extractfield(BootStatsIS_Q,'beta'))';                
                BootStatsOOS_Q = BootOOS(XBootOOS_Q, YBoot_Q, OOSBeg(2), PredSgn, gam, MeanVar(2), MeanRf(2));

                % run annual stuff - Have to bootstrap data again as the X variable is different
                XIS = VARS(vari).DataIS{3}; X = XIS(TBeg(1):TEnd(1)); Rm = RM(TBeg(1):TEnd(1)); Rf = RF(TBeg(1):TEnd(1));
                [XBootIS_M, RmBoot_M, RfBoot_M] = BootGenDataAR1(X, Rm, Rf, NBoot);
                T = size(XBootIS_M,1); XBootOOS_M = reshape(repmat(XBootIS_M, [T 1]), T, T, NBoot);
                XBootIS_ADec = XBootIS_M(12:12:end, :);
                XBootOOS_ADec = XBootOOS_M(12:12:end,12:12:end,:);
                T = size(XBootIS_ADec,1); RmBoot_ADec = NaN(T, NBoot); RfBoot_ADec = NaN(T, NBoot);
                for t = 1:T
                    RmBoot_ADec(t, :) = prod(1+RmBoot_M(12*t-11:12*t, :)) - 1;
                    RfBoot_ADec(t, :) = prod(1+RfBoot_M(12*t-11:12*t, :)) - 1;
                end
                YBoot_ADec = RmBoot_ADec - RfBoot_ADec;
                BootStatsIS_ADec = BootIS(XBootIS_ADec, YBoot_ADec);
                PredSgn = sign(extractfield(BootStatsIS_ADec,'beta'))';
                BootStatsOOS_ADec = BootOOS(XBootOOS_ADec, YBoot_ADec, OOSBeg(3), PredSgn, gam, MeanVar(3), MeanRf(3));

                % run annualJun stuff (using same bootstrapped annual data)
                XBootIS_AJun = XBootIS_ADec(1:end-1, :);
                XBootOOS_AJun = XBootOOS_ADec(1:end-1, 1:end-1, :);
                T = size(XBootIS_AJun,1); RmBoot_AJun = NaN(T, NBoot); RfBoot_AJun = NaN(T, NBoot);
                for t = 1:T
                    RmBoot_AJun(t, :) = prod(1+RmBoot_M(12*t-5:12*t+6, :)) - 1;
                    RfBoot_AJun(t, :) = prod(1+RfBoot_M(12*t-5:12*t+6, :)) - 1;
                end
                YBoot_AJun = RmBoot_AJun - RfBoot_AJun;
                BootStatsIS_AJun = BootIS(XBootIS_AJun, YBoot_AJun);
                PredSgn = sign(extractfield(BootStatsIS_AJun,'beta'))';                
                BootStatsOOS_AJun = BootOOS(XBootOOS_AJun, YBoot_AJun, OOSBeg(4), PredSgn, gam, MeanVar(4), MeanRf(4));

            elseif vari==7 % seminnual "skew"
                BootStatsIS_M = repmat(BootISstruct,NBoot,1); BootStatsOOS_M = repmat(BootOOSstruct,NBoot,1);
                BootStatsIS_Q = repmat(BootISstruct,NBoot,1); BootStatsOOS_Q = repmat(BootOOSstruct,NBoot,1);

                [TBeg, TEnd, OOSBeg, MeanVar, MeanRf] = CreateParams(vari,'AnnualDec');
                RM = RMarket_ADec; RF = RFree_ADec;
                XIS = VARS(vari).DataIS(2:2:end); X = XIS(TBeg(1):TEnd(1)); Rm = RM(TBeg(1):TEnd(1)); Rf = RF(TBeg(1):TEnd(1));
                [XBootIS_ADec, RmBoot_ADec, RfBoot_ADec] = BootGenDataAR1(X, Rm, Rf, NBoot);
                T = size(XBootIS_ADec,1); XBootOOS_ADec = reshape(repmat(XBootIS_ADec, [T 1]), T, T, NBoot);
                [BootStatsIS_ADec, BootStatsOOS_ADec] = RunAnnual(XBootIS_ADec, XBootOOS_ADec, RmBoot_ADec, RfBoot_ADec, OOSBeg, gam, MeanVar, MeanRf);

                % Need to do it all over again for June
                [TBeg, TEnd, OOSBeg, MeanVar, MeanRf] = CreateParams(vari,'AnnualJun');
                RM = RMarket_AJun; RF = RFree_AJun;
                XIS = VARS(vari).DataIS(2:2:end); X = XIS(TBeg(1)+1:TEnd(1)); Rm = RM(TBeg(1):TEnd(1)-1); Rf = RF(TBeg(1):TEnd(1)-1); % check leads/lags
                [XBootIS_AJun, RmBoot_AJun, RfBoot_AJun] = BootGenDataAR1(X, Rm, Rf, NBoot);
                T = size(XBootIS_AJun, 1); XBootOOS_AJun = reshape(repmat(XBootIS_AJun, [T 1]), T, T, NBoot);
                [BootStatsIS_AJun, BootStatsOOS_AJun] = RunAnnual(XBootIS_AJun, XBootOOS_AJun, RmBoot_AJun, RfBoot_AJun, OOSBeg, gam, MeanVar, MeanRf);

            elseif strcmp(VARS(vari).Freq,'Monthly')
                [TBeg, TEnd, OOSBeg, MeanVar, MeanRf] = CreateParams(vari,'Monthly');
                RM = RMarket_M; RF = RFree_M;
                XIS = VARS(vari).DataIS; X = XIS(TBeg(1):TEnd(1)); Rm = RM(TBeg(1):TEnd(1)); Rf = RF(TBeg(1):TEnd(1));
                [XBootIS_M, RmBoot_M, RfBoot_M] = BootGenDataAR1(X, Rm, Rf, NBoot);
                T = size(XBootIS_M,1); XBootOOS_M = reshape(repmat(XBootIS_M, [T 1]), T, T, NBoot);
                [BootStatsIS_M, BootStatsOOS_M, BootStatsIS_Q, BootStatsOOS_Q, BootStatsIS_ADec, BootStatsOOS_ADec, BootStatsIS_AJun, BootStatsOOS_AJun] = ...
                    RunMonthly(XBootIS_M, XBootOOS_M, RmBoot_M, RfBoot_M, OOSBeg, gam, MeanVar, MeanRf, NBoot);

            elseif strcmp(VARS(vari).Freq,'Quarterly')
                BootStatsIS_M = repmat(BootISstruct,NBoot,1); BootStatsOOS_M = repmat(BootOOSstruct,NBoot,1);

                [TBeg, TEnd, OOSBeg, MeanVar, MeanRf] = CreateParams(vari,'Quarterly');
                RM = RMarket_Q; RF = RFree_Q;
                XIS = VARS(vari).DataIS; X = XIS(TBeg(1):TEnd(1)); Rm = RM(TBeg(1):TEnd(1)); Rf = RF(TBeg(1):TEnd(1));
                [XBootIS_Q, RmBoot_Q, RfBoot_Q] = BootGenDataAR1(X, Rm, Rf, NBoot);
                T = size(XBootIS_Q,1); XBootOOS_Q = reshape(repmat(XBootIS_Q, [T 1]), T, T, NBoot);
                [BootStatsIS_Q, BootStatsOOS_Q, BootStatsIS_ADec, BootStatsOOS_ADec, BootStatsIS_AJun, BootStatsOOS_AJun] = ...
                    RunQuarterly(XBootIS_Q, XBootOOS_Q, RmBoot_Q, RfBoot_Q, OOSBeg, gam, MeanVar, MeanRf, NBoot);

            elseif strcmp(VARS(vari).Freq,'Annual') && (vari==23 || vari==24) % DGP is NOT AR1
                BootStatsIS_M = repmat(BootISstruct,NBoot,1); BootStatsOOS_M = repmat(BootOOSstruct,NBoot,1);
                BootStatsIS_Q = repmat(BootISstruct,NBoot,1); BootStatsOOS_Q = repmat(BootOOSstruct,NBoot,1);

                [TBeg, TEnd, OOSBeg, MeanVar, MeanRf] = CreateParams(vari,'AnnualDec');
                RM = RMarket_ADec; RF = RFree_ADec;
                XIS = VARS(vari).DataIS; X = XIS(TBeg(1):TEnd(1)); Rm = RM(TBeg(1):TEnd(1)); Rf = RF(TBeg(1):TEnd(1));
                [XBootIS_ADec, RmBoot_ADec, RfBoot_ADec] = BootGenData(X, Rm, Rf, NBoot);
                T = size(XBootIS_ADec,1); XBootOOS_ADec = reshape(repmat(XBootIS_ADec, [T 1]), T, T, NBoot);
                [BootStatsIS_ADec, BootStatsOOS_ADec] = RunAnnual(XBootIS_ADec, XBootOOS_ADec, RmBoot_ADec, RfBoot_ADec, OOSBeg, gam, MeanVar, MeanRf);

                % Need to do it all over again for June
                [TBeg, TEnd, OOSBeg, MeanVar, MeanRf] = CreateParams(vari,'AnnualJun');
                RM = RMarket_AJun; RF = RFree_AJun;
                XIS = VARS(vari).DataIS; X = XIS(TBeg(1)+1:TEnd(1)); Rm = RM(TBeg(1):TEnd(1)-1); Rf = RF(TBeg(1):TEnd(1)-1); % check leads/lags
                [XBootIS_AJun, RmBoot_AJun, RfBoot_AJun] = BootGenData(X, Rm, Rf, NBoot);
                T = size(XBootIS_AJun, 1); XBootOOS_AJun = reshape(repmat(XBootIS_AJun, [T 1]), T, T, NBoot);
                [BootStatsIS_AJun, BootStatsOOS_AJun] = RunAnnual(XBootIS_AJun, XBootOOS_AJun, RmBoot_AJun, RfBoot_AJun, OOSBeg, gam, MeanVar, MeanRf);

            elseif strcmp(VARS(vari).Freq,'Annual')
                BootStatsIS_M = repmat(BootISstruct,NBoot,1); BootStatsOOS_M = repmat(BootOOSstruct,NBoot,1);
                BootStatsIS_Q = repmat(BootISstruct,NBoot,1); BootStatsOOS_Q = repmat(BootOOSstruct,NBoot,1);

                [TBeg, TEnd, OOSBeg, MeanVar, MeanRf] = CreateParams(vari,'AnnualDec');
                RM = RMarket_ADec; RF = RFree_ADec;
                XIS = VARS(vari).DataIS; X = XIS(TBeg(1):TEnd(1)); Rm = RM(TBeg(1):TEnd(1)); Rf = RF(TBeg(1):TEnd(1));
                [XBootIS_ADec, RmBoot_ADec, RfBoot_ADec] = BootGenDataAR1(X, Rm, Rf, NBoot);
                T = size(XBootIS_ADec,1); XBootOOS_ADec = reshape(repmat(XBootIS_ADec, [T 1]), T, T, NBoot);
                [BootStatsIS_ADec, BootStatsOOS_ADec] = RunAnnual(XBootIS_ADec, XBootOOS_ADec, RmBoot_ADec, RfBoot_ADec, OOSBeg, gam, MeanVar, MeanRf);

                % Need to do it all over again for June
                [TBeg, TEnd, OOSBeg, MeanVar, MeanRf] = CreateParams(vari,'AnnualJun');
                RM = RMarket_AJun; RF = RFree_AJun;
                XIS = VARS(vari).DataIS; X = XIS(TBeg(1)+1:TEnd(1)); Rm = RM(TBeg(1):TEnd(1)-1); Rf = RF(TBeg(1):TEnd(1)-1); % check leads/lags
                [XBootIS_AJun, RmBoot_AJun, RfBoot_AJun] = BootGenDataAR1(X, Rm, Rf, NBoot);
                T = size(XBootIS_AJun, 1); XBootOOS_AJun = reshape(repmat(XBootIS_AJun, [T 1]), T, T, NBoot);
                [BootStatsIS_AJun, BootStatsOOS_AJun] = RunAnnual(XBootIS_AJun, XBootOOS_AJun, RmBoot_AJun, RfBoot_AJun, OOSBeg, gam, MeanVar, MeanRf);
            end

            AllBootStatsIS_M(:,vari) = BootStatsIS_M;   AllBootStatsIS_Q (:,vari) = BootStatsIS_Q;  AllBootStatsIS_ADec(:,vari) = BootStatsIS_ADec;   AllBootStatsIS_AJun(:,vari) = BootStatsIS_AJun;
            AllBootStatsOOS_M(:,vari) = BootStatsOOS_M; AllBootStatsOOS_Q(:,vari) = BootStatsOOS_Q; AllBootStatsOOS_ADec(:,vari) = BootStatsOOS_ADec; AllBootStatsOOS_AJun(:,vari) = BootStatsOOS_AJun;
        end

        save(strcat('BootAllResults',num2str(BootSample),num2str(SimpleLog)),'AllBoot*');
    end % end for SimpleLog = 1:2
end % end for BootSample = 1:3
return


%%
function [TBeg, TEnd, OOSBeg, MeanVar, MeanRf] = CreateParams(vari,Type)
global RMarket_M RMarket_Q RMarket_ADec RMarket_AJun RFree_M RFree_Q RFree_ADec RFree_AJun MktVar_M MktVar_Q MktVar_ADec MktVar_AJun yyyymm yyyyq yyyy SimpleLog BootSample BootISstruct BootOOSstruct rootdir;

load(strcat(rootdir,'Data\P00Data'));
if BootSample==1
    load(strcat(rootdir,'Analysis\Sample1_Full\ResultsFile'));
elseif BootSample==3
    load(strcat(rootdir,'Analysis\Sample3_BegAuthor_EndFull\ResultsFile'));
end

if BootSample==2 %Only AuthorSample
    for n = 1:46
        VARS(n).SampleBeg_M = VARS(n).SampleBegAuthor_M; VARS(n).SampleEnd_M = VARS(n).SampleEndAuthor_M; VARS(n).OOSBeg_M = VARS(n).OOSBegAuthor_M;
        VARS(n).SampleBeg_Q = VARS(n).SampleBegAuthor_Q; VARS(n).SampleEnd_Q = VARS(n).SampleEndAuthor_Q; VARS(n).OOSBeg_Q = VARS(n).OOSBegAuthor_Q;
        VARS(n).SampleBeg_A = VARS(n).SampleBegAuthor_A; VARS(n).SampleEnd_A = VARS(n).SampleEndAuthor_A; VARS(n).OOSBeg_A = VARS(n).OOSBegAuthor_A;
    end
elseif BootSample==3 % From AuthorBeg to SampleEnd
    for n = 1:46
        VARS(n).SampleBeg_M = VARS(n).SampleBegAuthor_M; VARS(n).OOSBeg_M = VARS(n).OOSBegAuthor_M;
        VARS(n).SampleBeg_Q = VARS(n).SampleBegAuthor_Q; VARS(n).OOSBeg_Q = VARS(n).OOSBegAuthor_Q;
        VARS(n).SampleBeg_A = VARS(n).SampleBegAuthor_A; VARS(n).OOSBeg_A = VARS(n).OOSBegAuthor_A;
    end
end

if strcmp(Type,'Monthly')
    TBeg(1) = find(yyyymm==VARS(vari).SampleBeg_M); TEnd(1) = find(yyyymm==VARS(vari).SampleEnd_M); OOSBeg(1) = find(yyyymm==VARS(vari).OOSBeg_M);
    TBeg(2) = find(yyyyq==VARS(vari).SampleBeg_Q); TEnd(2) = find(yyyyq==VARS(vari).SampleEnd_Q); OOSBeg(2) = find(yyyyq==VARS(vari).OOSBeg_Q);
    TBeg(3) = find(yyyy==VARS(vari).SampleBeg_A); TEnd(3) = find(yyyy==VARS(vari).SampleEnd_A); OOSBeg(3) = find(yyyy==VARS(vari).OOSBeg_A);
    TBeg(4) = TBeg(3) + 1; TEnd(4) = TEnd(3); OOSBeg(4) = OOSBeg(3) + 1;

    MeanVar = [mean(MktVar_M(OOSBeg(1)-1:TEnd(1)-1)); mean(MktVar_Q(OOSBeg(2)-1:TEnd(2)-1)); mean(MktVar_ADec(OOSBeg(3)-1:TEnd(3)-1)); mean(MktVar_AJun(OOSBeg(4)-1:TEnd(4)-1))];
    MeanRf = [mean(RFree_M(OOSBeg(1):TEnd(1))); mean(RFree_Q(OOSBeg(2):TEnd(2))); mean(RFree_ADec(OOSBeg(3):TEnd(3))); mean(RFree_AJun(OOSBeg(4):TEnd(4)))];
    OOSBeg = OOSBeg - TBeg + 1;
elseif strcmp(Type,'Quarterly')
    TBeg(1) = find(yyyyq==VARS(vari).SampleBeg_Q); TEnd(1) = find(yyyyq==VARS(vari).SampleEnd_Q); OOSBeg(1) = find(yyyyq==VARS(vari).OOSBeg_Q);
    TBeg(2) = find(yyyy==VARS(vari).SampleBeg_A); TEnd(2) = find(yyyy==VARS(vari).SampleEnd_A); OOSBeg(2) = find(yyyy==VARS(vari).OOSBeg_A);
    TBeg(3) = TBeg(2) + 1; TEnd(3) = TEnd(2); OOSBeg(3) = OOSBeg(2) + 1;

    MeanVar = [mean(MktVar_Q(OOSBeg(1)-1:TEnd(1)-1)); mean(MktVar_ADec(OOSBeg(2)-1:TEnd(2)-1)); mean(MktVar_AJun(OOSBeg(3)-1:TEnd(3)-1))];
    MeanRf = [mean(RFree_Q(OOSBeg(1):TEnd(1))); mean(RFree_ADec(OOSBeg(2):TEnd(2))); mean(RFree_AJun(OOSBeg(3):TEnd(3)))];
    OOSBeg = OOSBeg - TBeg + 1;
elseif strcmp(Type,'AnnualDec')
    TBeg = find(yyyy==VARS(vari).SampleBeg_A); TEnd = find(yyyy==VARS(vari).SampleEnd_A); OOSBeg = find(yyyy==VARS(vari).OOSBeg_A);

    MeanVar = mean(MktVar_ADec(OOSBeg-1:TEnd-1));
    MeanRf = mean(RFree_ADec(OOSBeg:TEnd));
    OOSBeg = OOSBeg - TBeg + 1;
elseif strcmp(Type,'AnnualJun')
    TBeg = find(yyyy==VARS(vari).SampleBeg_A); TEnd = find(yyyy==VARS(vari).SampleEnd_A); OOSBeg = find(yyyy==VARS(vari).OOSBeg_A);
    TBeg = TBeg + 1; OOSBeg = OOSBeg + 1;

    MeanVar = mean(MktVar_AJun(OOSBeg-1:TEnd-1));
    MeanRf = mean(RFree_AJun(OOSBeg:TEnd));
    OOSBeg = OOSBeg - TBeg + 1;
end
return

%%
function [BootStatsIS_M, BootStatsOOS_M, BootStatsIS_Q, BootStatsOOS_Q, BootStatsIS_ADec, BootStatsOOS_ADec, BootStatsIS_AJun, BootStatsOOS_AJun] = ...
    RunMonthly(XBootIS_M, XBootOOS_M, RmBoot_M, RfBoot_M, OOSBeg, gam, MeanVar, MeanRf, NBoot)

% run monthly stuff
YBoot_M = RmBoot_M - RfBoot_M;
BootStatsIS_M = BootIS(XBootIS_M, YBoot_M);
PredSgn = sign(extractfield(BootStatsIS_M,'beta'))';
BootStatsOOS_M = BootOOS(XBootOOS_M, YBoot_M, OOSBeg(1), PredSgn, gam, MeanVar(1), MeanRf(1));

% run quarterly stuff (using same bootstrapped data)
XBootIS_Q = XBootIS_M(3:3:end, :);
XBootOOS_Q = XBootOOS_M(3:3:end,3:3:end,:);
T = size(XBootIS_Q,1); RmBoot_Q = NaN(T, NBoot); RfBoot_Q = NaN(T, NBoot);
for t = 1:T
    RmBoot_Q(t, :) = prod(1+RmBoot_M(3*t-2:3*t, :)) - 1;
    RfBoot_Q(t, :) = prod(1+RfBoot_M(3*t-2:3*t, :)) - 1;
end
YBoot_Q = RmBoot_Q - RfBoot_Q;
BootStatsIS_Q = BootIS(XBootIS_Q, YBoot_Q);
PredSgn = sign(extractfield(BootStatsIS_Q,'beta'))';
BootStatsOOS_Q = BootOOS(XBootOOS_Q, YBoot_Q, OOSBeg(2), PredSgn, gam, MeanVar(2), MeanRf(2));

% run annualDec stuff (using same bootstrapped data)
XBootIS_ADec = XBootIS_M(12:12:end, :);
XBootOOS_ADec = XBootOOS_M(12:12:end,12:12:end,:);
T = size(XBootIS_ADec,1); RmBoot_ADec = NaN(T, NBoot); RfBoot_ADec = NaN(T, NBoot);
for t = 1:T
    RmBoot_ADec(t, :) = prod(1+RmBoot_M(12*t-11:12*t, :)) - 1;
    RfBoot_ADec(t, :) = prod(1+RfBoot_M(12*t-11:12*t, :)) - 1;
end
% Data are aligned so that XBootIS_M(12) lines up with RBoot_M(1:12), etc.
YBoot_ADec = RmBoot_ADec - RfBoot_ADec;
BootStatsIS_ADec = BootIS(XBootIS_ADec, YBoot_ADec);
PredSgn = sign(extractfield(BootStatsIS_ADec,'beta'))';
BootStatsOOS_ADec = BootOOS(XBootOOS_ADec, YBoot_ADec, OOSBeg(3), PredSgn, gam, MeanVar(3), MeanRf(3));

% run annualJun stuff (using same bootstrapped data)
XBootIS_AJun = XBootIS_ADec(1:end-1, :);
XBootOOS_AJun = XBootOOS_ADec(1:end-1, 1:end-1, :);
T = size(XBootIS_AJun,1); RmBoot_AJun = NaN(T, NBoot); RfBoot_AJun = NaN(T, NBoot);
for t = 1:T
    RmBoot_AJun(t, :) = prod(1+RmBoot_M(12*t-5:12*t+6, :)) - 1;
    RfBoot_AJun(t, :) = prod(1+RfBoot_M(12*t-5:12*t+6, :)) - 1;
end
% Data are aligned so that XBootIS_M(12) lines up with RBoot_M(7:18), etc.
YBoot_AJun = RmBoot_AJun - RfBoot_AJun;
BootStatsIS_AJun = BootIS(XBootIS_AJun, YBoot_AJun);
PredSgn = sign(extractfield(BootStatsIS_AJun,'beta'))';
BootStatsOOS_AJun = BootOOS(XBootOOS_AJun, YBoot_AJun, OOSBeg(4), PredSgn, gam, MeanVar(4), MeanRf(4));
return

%%
function [BootStatsIS_Q, BootStatsOOS_Q, BootStatsIS_ADec, BootStatsOOS_ADec, BootStatsIS_AJun, BootStatsOOS_AJun] = ...
    RunQuarterly(XBootIS_Q, XBootOOS_Q, RmBoot_Q, RfBoot_Q, OOSBeg, gam, MeanVar, MeanRf, NBoot)

% run quarterly stuff
YBoot_Q = RmBoot_Q - RfBoot_Q;
BootStatsIS_Q = BootIS(XBootIS_Q, YBoot_Q);
PredSgn = sign(extractfield(BootStatsIS_Q,'beta'))';
BootStatsOOS_Q = BootOOS(XBootOOS_Q, YBoot_Q, OOSBeg(1), PredSgn, gam, MeanVar(1), MeanRf(1));

% run annualDec stuff (using same bootstrapped data)
XBootIS_ADec = XBootIS_Q(4:4:end, :);
XBootOOS_ADec = XBootOOS_Q(4:4:end,4:4:end,:);
T = size(XBootIS_ADec,1); RmBoot_ADec = NaN(T, NBoot); RfBoot_ADec = NaN(T, NBoot);
for t = 1:T
    RmBoot_ADec(t, :) = prod(1+RmBoot_Q(4*t-3:4*t, :)) - 1;
    RfBoot_ADec(t, :) = prod(1+RfBoot_Q(4*t-3:4*t, :)) - 1;
end
YBoot_ADec = RmBoot_ADec - RfBoot_ADec;
BootStatsIS_ADec = BootIS(XBootIS_ADec, YBoot_ADec);
PredSgn = sign(extractfield(BootStatsIS_ADec,'beta'))';
BootStatsOOS_ADec = BootOOS(XBootOOS_ADec, YBoot_ADec, OOSBeg(2), PredSgn, gam, MeanVar(2), MeanRf(2));

% run annualJun stuff (using same bootstrapped data)
XBootIS_AJun = XBootIS_ADec(1:end-1, :);
XBootOOS_AJun = XBootOOS_ADec(1:end-1, 1:end-1, :);
T = size(XBootIS_AJun, 1); RmBoot_AJun = NaN(T, NBoot); RfBoot_AJun = NaN(T, NBoot);
for t = 1:T
    RmBoot_AJun(t, :) = prod(1+RmBoot_Q(4*t-1:4*t+2, :)) - 1;
    RfBoot_AJun(t, :) = prod(1+RfBoot_Q(4*t-1:4*t+2, :)) - 1;
end
YBoot_AJun = RmBoot_AJun - RfBoot_AJun;
BootStatsIS_AJun = BootIS(XBootIS_AJun, YBoot_AJun);
PredSgn = sign(extractfield(BootStatsIS_AJun,'beta'))';
BootStatsOOS_AJun = BootOOS(XBootOOS_AJun, YBoot_AJun, OOSBeg(3), PredSgn, gam, MeanVar(3), MeanRf(3));
return

%%
function [BootStatsIS_A, BootStatsOOS_A] = RunAnnual(XBootIS_A, XBootOOS_A, RmBoot_A, RfBoot_A, OOSBeg, gam, MeanVar, MeanRf)

YBoot_A = RmBoot_A - RfBoot_A;
BootStatsIS_A = BootIS(XBootIS_A, YBoot_A);
PredSgn = sign(extractfield(BootStatsIS_A,'beta'))';
BootStatsOOS_A = BootOOS(XBootOOS_A, YBoot_A, OOSBeg, PredSgn, gam, MeanVar, MeanRf);
return

%%
function [BootStatsIS_M, BootStatsOOS_M, BootStatsIS_Q, BootStatsOOS_Q, BootStatsIS_ADec, BootStatsOOS_ADec, BootStatsIS_AJun, BootStatsOOS_AJun] = ...
    RunMonthlyPLS(XPrePLSBoot_M, RmBoot_M, RfBoot_M, OOSBeg, gam, MeanVar, MeanRf, NBoot)

T = size(RmBoot_M, 1);

% run monthly stuff
YBoot_M = RmBoot_M - RfBoot_M;
XBootIS_M = NaN(T, NBoot); XBootOOS_M = NaN(T, T, NBoot);
for b = 1:NBoot
    XBootOOS_M(:,:,b) = PLS(squeeze(XPrePLSBoot_M(:,b,:)),  YBoot_M(:,b), 1, 120, T);
    XBootIS_M(:,b) = XBootOOS_M(:,end,b);
end
BootStatsIS_M = BootIS(XBootIS_M, YBoot_M);
PredSgn = sign(extractfield(BootStatsIS_M,'beta'))';
BootStatsOOS_M = BootOOS(XBootOOS_M, YBoot_M, OOSBeg(1), PredSgn, gam, MeanVar(1), MeanRf(1));

% run quarterly stuff (using same bootstrapped XPrePLSBoot_M but redo PLS)
TQ = floor(T/3);
RmBoot_Q = NaN(TQ, NBoot); RfBoot_Q = NaN(TQ, NBoot);
for t = 1:TQ
    RmBoot_Q(t, :) = prod(1+RmBoot_M(3*t-2:3*t, :)) - 1;
    RfBoot_Q(t, :) = prod(1+RfBoot_M(3*t-2:3*t, :)) - 1;
end
YBoot_Q = RmBoot_Q - RfBoot_Q;
XBootIS_Q = NaN(TQ, NBoot); XBootOOS_Q = NaN(TQ, TQ, NBoot);
for b = 1:NBoot
    XBootOOS_Q(:,:,b) = PLS(squeeze(XPrePLSBoot_M(3:3:end,b,:)), YBoot_Q(:,b), 1, 40, TQ);
    XBootIS_Q(:,b) = XBootOOS_Q(:,end,b);
end
BootStatsIS_Q = BootIS(XBootIS_Q, YBoot_Q);
PredSgn = sign(extractfield(BootStatsIS_Q,'beta'))';
BootStatsOOS_Q = BootOOS(XBootOOS_Q, YBoot_Q, OOSBeg(2), PredSgn, gam, MeanVar(2), MeanRf(2));

% run annualDec stuff (using same bootstrapped XPrePLSBoot_M but redo PLS)
TA = floor(T/12);
RmBoot_ADec = NaN(TA, NBoot); RfBoot_ADec = NaN(TA, NBoot);
for t = 1:TA
    RmBoot_ADec(t, :) = prod(1+RmBoot_M(12*t-11:12*t, :)) - 1;
    RfBoot_ADec(t, :) = prod(1+RfBoot_M(12*t-11:12*t, :)) - 1;
end
YBoot_ADec = RmBoot_ADec - RfBoot_ADec;
XBootIS_ADec = NaN(TA, NBoot); XBootOOS_ADec = NaN(TA, TA, NBoot);
for b = 1:NBoot
    XBootOOS_ADec(:,:,b) = PLS(squeeze(XPrePLSBoot_M(12:12:end,b,:)), YBoot_ADec(:,b), 1, 10, TA);
    XBootIS_ADec(:,b) = XBootOOS_ADec(:,end,b);
end
% Data are aligned so that XPreIS_M(12) lines up with YBoot_M(1:12), etc.
BootStatsIS_ADec = BootIS(XBootIS_ADec, YBoot_ADec);
PredSgn = sign(extractfield(BootStatsIS_ADec,'beta'))';
BootStatsOOS_ADec = BootOOS(XBootOOS_ADec, YBoot_ADec, OOSBeg(3), PredSgn, gam, MeanVar(3), MeanRf(3));

% run annualJun stuff (using same bootstrapped XPrePLSBoot_M but redo PLS)
TA = floor(T/12);
if mod(T,12)==0
    TA = TA - 1;
    XPrePLSBoot_M(end-11:end,:,:) = [];
end
RmBoot_AJun = NaN(TA, NBoot); RfBoot_AJun = NaN(TA, NBoot);
for t = 1:TA
    RmBoot_AJun(t, :) = prod(1+RmBoot_M(12*t-5:12*t+6, :)) - 1;
    RfBoot_AJun(t, :) = prod(1+RfBoot_M(12*t-5:12*t+6, :)) - 1;
end
YBoot_AJun = RmBoot_AJun - RfBoot_AJun;
XBootIS_AJun = NaN(TA, NBoot); XBootOOS_AJun = NaN(TA, TA, NBoot);
for b = 1:NBoot
    XBootOOS_AJun(:,:,b) = PLS(squeeze(XPrePLSBoot_M(12:12:end,b,:)), YBoot_AJun(:,b), 1, 10, TA);
    XBootIS_AJun(:,b) = XBootOOS_AJun(:,end,b);
end
% Data are aligned so that XPreIS_M(12) lines up with YBoot_M(7:18), etc.
BootStatsIS_AJun = BootIS(XBootIS_AJun, YBoot_AJun);
PredSgn = sign(extractfield(BootStatsIS_AJun,'beta'))';
BootStatsOOS_AJun = BootOOS(XBootOOS_AJun, YBoot_AJun, OOSBeg(4), PredSgn, gam, MeanVar(4), MeanRf(4));
return
